#!/bin/bash
set -e

while getopts  ":s:p:" opts
do
        case  $opts  in
        s)
		barcode=$OPTARG
		sample_name=$OPTARG
		
		;;
		p)
		out_prefix=$OPTARG
		
		;;
		?)
		echo `basename $0` "<input.vcf> <out.txt>"
		exit 1
		;;
        esac
done
shift $(($OPTIND - 1))

if [ -z "$2" ]; then
	echo `basename $0` "<input.vcf> <out.txt>"
	exit 1
fi


#-----------------------------------------------
. /mnt/ilustre/app/medical/tools/.var #---------
#-----------------------------------------------


dat=`date +'%b_%d_%H_%M_%S_%Y'`

fil=`basename $1`

dir_name=${dat}__$fil

mkdir $dir_name
echo $dir_name maked 
echo tep files will be located in $dir_name ...
cd $dir_name


script_path=$tools_path/script
snpeff_path=${tools_path}/snpEff
snpsift=${snpeff_path}/SnpSift.jar

java_memory=1g
echo java memory: $java_memory



echo
echo extractField 
echo

echo CHROM
java -jar $snpsift extractField ../$1 CHROM > chrom.txt

echo POS
java -jar $snpsift extractField ../$1 POS > pos.txt

# 将 ID 域预处理，然后将其分离为单独的列
java -jar $snpsift extractField ../$1 ID | sed -e 's/^COSM/;COSM/' | sed -e 's/\(^rs[0-9]*$\)/\1;/' > 1.id.sep.txt

echo RS ID
cut -d ";" -f 1 1.id.sep.txt > rs_id.txt
echo COSMIC ID
cut -d ";" -f 2 1.id.sep.txt > cosmic_id.txt

echo REF
java -Xmx$java_memory -jar $snpsift extractField ../$1 REF > 1.ref.txt
echo ALT
java -Xmx$java_memory -jar $snpsift extractField ../$1 ALT > 1.alt.txt
cp 1.ref.txt ref.txt
cp 1.alt.txt alt.txt



# DP REF ALT 需要计算或需要再格式化，所以都要单独提取
# java -Xmx$java_memory -jar $snpsift extractField ../$1 "GEN[0].DP" "DP" | sed -e 's/\t/\//' > 1.dp2.txt
echo DP
java -Xmx$java_memory -jar $snpsift extractField ../$1 "GEN[0].DP" > 1.dp.txt
cp 1.dp.txt dp.txt




# 取得 ref alt 同位基因的深度值 such as 23,22,12
# then merge like this ref=#dp;alt=#dp
echo
echo extractField ../$1 GEN[0].AD[*]
java -Xmx$java_memory -jar $snpsift extractField ../$1 "GEN[0].AD[*]" > 1.ref.alt.dp.txt
paste 1.ref.alt.dp.txt 1.dp.txt > 1.all_dp.txt
java -Xmx$java_memory -jar $snpsift extractField ../$1 "GEN[0].AD[0]" > 1.ref.dp.txt
# -s 保存文件时的域分隔符 -d 从第二位开始取，即为 所有的alt
java -Xmx$java_memory -jar $snpsift extractField ../$1 -s "," "GEN[0].AD[*]" | cut -d "," -f 2- > 1.alt.dp.txt

# format like this: ref=#dp alt1,alt2...,atln=#dp1,#dp2...,#dpn
# then merge like this: ref;alt
paste -d "=" 1.ref.txt  1.ref.dp.txt > 1.ref.dp.merge.txt
paste -d "=" 1.alt.txt  1.alt.dp.txt > 1.alt.dp.merge.txt
paste -d ";" 1.ref.dp.merge.txt 1.alt.dp.merge.txt > 1.all.dp.txt

echo ALLELE COVERAGE
cp 1.all.dp.txt allele_cvg.txt




# 删除第一行的注释用于计算
sed -i "1d" 1.all_dp.txt
perl $script_path/format_0318.pl 1.all_dp.txt > 1.all.frq.txt



# 取得 ref 的频率并 format like this: ref=#ref_frq
cut -f 1 1.all.frq.txt > 1.ref.frq.txt
paste -d "=" 1.ref.txt 1.ref.frq.txt > 1.ref_frq.txt

# 取得 alt 的频率并 format like this: alt1,alt2=#alt1_frq,alt2_frq
cut -f 2 1.all.frq.txt > 1.alt1.frq.txt




cut -f 2- 1.all.frq.txt | sed -e 's/\t/,/' > 1.alt.frq.txt
paste -d "=" 1.alt.txt 1.alt.frq.txt > 1.alt_frq.txt
paste -d ";" 1.ref_frq.txt 1.alt_frq.txt > 1.all_frq.txt
echo ALT FREQUENCY
# 在第一行插入注释 用 echo cat 而不用 sed 1i 或 1c 是因为如果 file 为空，可能不起作用
echo ALT FREQUENCY |cat - 1.alt1.frq.txt > alt_freq.txt

# 计算 p-value ??
echo
echo QUAL
java -Xmx$java_memory -jar $snpsift extractField ../$1 QUAL > 1.qual.txt
sed -i "1d" 1.qual.txt
echo
echo format p-value
perl $script_path/format_0318.pl 1.qual.txt > 1.p.txt
# sed -i "1i QUAL" 1.qual.txt
# sed -i "1i P_value" 1.p.txt

echo P-VALUE
echo P-VALUE |cat - 1.p.txt > p_value.txt


echo
echo VARTYPE
java -Xmx$java_memory -jar $snpsift extractField ../$1 VARTYPE > vartype.txt

# cut -f1-2 i.txt >i1.txt
# cut -f3 i.txt | enst_refseq.pl -|awk '{print $1 "/" $2}' > i2_trid.txt
# cut -f4- i.txt > i3.txt
# paste i1.txt i2_trid.txt i3.txt > i.new.txt



# paste 4.txt 1.ref.txt 1.alt.txt i.new.txt 1.dp2.txt 1.all.dp.txt 1.all_frq.txt 1.alt1.frq.txt 1.p.txt > ii.txt



echo FILTER
java -Xmx$java_memory -jar $snpsift extractField ../$1 FILTER > filter.txt

echo INFO
format_vep.sh ../$1 1.csq.txt
format_vep_info.pl 1.csq.txt > hgvs.txt

cat hgvs.txt |sed 's/<br \/>/    /g'> hgvs.xls.txt

# cat 1.csq.txt|cut -d'|' -f2 > consequence.txt
# cat 1.csq.txt|cut -d'|' -f3 > impact.txt
# cat 1.csq.txt|cut -d'|' -f4 > gene.txt
# cat 1.csq.txt|cut -d'|' -f7 > refseq.txt
# cat 1.csq.txt|cut -d'|' -f8 > biotype.txt
# cat 1.csq.txt|cut -d'|' -f9,10 |sed -e 's/|//' > exon.intron.txt
# cat 1.csq.txt|cut -d'|' -f11 > hgvsc.txt
# cat 1.csq.txt|cut -d'|' -f12 > hgvsp.txt

# 找到了才替换并打印，找不到的不打印
# cat 1.csq.txt |sed -n 's/,/<br \/>/pg'> info.txt
cat 1.csq.txt |sed 's/,/<br \/>/g'> info.txt
cat 1.csq.txt |sed 's/,/<br \/>/g'> info.html.txt
cat 1.csq.txt |sed 's/,/    /g'> info.xls.txt
# cat 1.csq.txt |sed 's/,/\\n/g'> info.xls.txt

echo Index
seq 1 $[`wc -l chrom.txt |cut -d' ' -f1`-1] >> index.txt

sed -i '1c\Chrom' chrom.txt
sed -i '1c\GRCh37 Pos' pos.txt
sed -i '1c\dbSNP ID' rs_id.txt
sed -i '1c\COSMIC ID' cosmic_id.txt
sed -i '1c\Ref' ref.txt
sed -i '1c\Var' alt.txt
sed -i '1c\VarType' vartype.txt
sed -i '1c\Coverage' dp.txt
sed -i '1c\Allele Coverage' allele_cvg.txt
sed -i '1c\Allele Frequency' alt_freq.txt
sed -i '1c\P-VALUE' p_value.txt
sed -i '1c\Filter' filter.txt
sed -i '1c\Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|DISTANCE|STRAND|SYMBOL_SOURCE|HGNC_ID|REFSEQ_MATCH|HGVS_OFFSET' info.txt
sed -i '1c\INFO' hgvs.txt
sed -i '1c\Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|DISTANCE|STRAND|SYMBOL_SOURCE|HGNC_ID|REFSEQ_MATCH|HGVS_OFFSET' info.html.txt

sed -i '1c\Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|DISTANCE|STRAND|SYMBOL_SOURCE|HGNC_ID|REFSEQ_MATCH|HGVS_OFFSET' info.xls.txt
sed -i '1c\INFO' hgvs.xls.txt



# sed -i '1c\Gene' gene.txt
# sed -i '1c\Transcript ID' refseq.txt
# sed -i '1c\Exon/Intron Rank' exon.intron.txt
# sed -i '1c\Consequence' consequence.txt
# sed -i '1c\HGVS_Codon Change' hgvsc.txt
# sed -i '1c\HGVS_AA Change' hgvsp.txt
# sed -i '1c\Impact' impact.txt



paste chrom.txt pos.txt rs_id.txt cosmic_id.txt ref.txt alt.txt vartype.txt dp.txt allele_cvg.txt alt_freq.txt p_value.txt filter.txt hgvs.xls.txt info.xls.txt> ../$2

paste chrom.txt pos.txt rs_id.txt cosmic_id.txt ref.txt alt.txt vartype.txt dp.txt allele_cvg.txt alt_freq.txt p_value.txt filter.txt hgvs.txt info.html.txt> ../$3

# for txt in chrom.txt pos.txt rs_id.txt cosmic_id.txt ref.txt alt.txt vartype.txt dp.txt allele_cvg.txt alt_freq.txt p_value.txt filter.txt hgvs.txt info.txt; do
# txt_base=`basename $txt .txt`
# cat $txt |awk '{print "\"" $0 "\""}' >$txt_base.csv
# done


# paste -d',' chrom.csv pos.csv rs_id.csv cosmic_id.csv ref.csv alt.csv vartype.csv dp.csv allele_cvg.csv alt_freq.csv p_value.csv filter.csv hgvs.csv info.csv> 1.csv

# cat 1.csv |awk '{gsub("<br />", "\012");print}' >../$3

# paste ii.txt iii.txt > iv.txt

# sed -i "1d" iv.txt
# cat $data_path/title1.txt iv.txt > ../$2



# for sift and polyphen2 format

# considering grchb37 and hg19 are the same in some extent...
# cat 2.txt |awk '{sub(/^/, "chr");sub(/\t/, ":"); print}' > 1.chr_pos.txt
# paste -d "/" 1.ref.txt 1.alt.txt > 1.ref_alt.txt
# paste -d " " 1.chr_pos.txt 1.ref_alt.txt | awk 'NR>1 {print}' > ../1.pph.txt










